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Abstract 

We discuss the general design of the ANTARES code which is intended for 
simulations in stellar hydrodynamics with radiative transfer and realistic micro- 
physics in ID, 2D and 3D. We then compare the quality of various numerical 
methods. We have applied ANTARES in order to obtain high resolution simu- 
lations of solar granulation which we describe and analyze. In order to obtain 
high resolution, we apply grid refinement to a region predominantly occupied by 
an exploding granule. Strong, rapidly rotating vortex tubes of small diameter 
(~ 100 km) generated by the downdrafts and ascending into the photosphere 
near the granule boundaries evolve, often entering the photosphere from below 
in an arclike fashion. They essentially contribute to the turbulent velocity field 
near the granule boundaries. 
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1. Introduction 

Modelling of hydrodynamic flows in the setting as appropriate for stellar 
physics has provided a large number of insights over the past decades. If we 
stick to the vicinityof solar granulation studies only, we may mention work such 

n, 0, [m, [si], [13]. 



as 



In order to conduct research in the context of stellar convection and related 
areas we are developing a code for such simulations. In this paper, we describe 



* corresponding author 
Email address: herbert .muthsamOunivie . ac . at (H.J. Muthsam) 



Preprint submitted to New Astronomy 



May 2, 2009 



the most central and most mature features of the continually expanding code, 
i.e. hydrodynamics with radiative transfer and realistic material properties for 
the stellar and of course solar case. We discuss the relative performance of 
various numerical methods. We exemplify the use and benefits granted by the 
program in discussing high resolution simulations of solar granulation. 

ANTARES (A Numerical Tool for Astrophysical RESearch) is designed ac- 
cording to the following general principles and aims: 

• general: time dependent compressible hydrodynamics and extensions 
(such as MHD) in ID, 2D and 3D; Fortran90 code with a modular struc- 
ture; 

• numerics: various high resolution numerical schemes of conservative form 
implemented in order to evaluate their merits; 

• radiative transfer: short characteristics method (use of the diffusion 
approximation where appropriate); either grey or nongrey by the binning 
method based on state of the art opacities; 

• microphysics: realistic (or idealized); 

• gridding: logically rectangular; either rectangular or spherical coordi- 
nates; equidistant or (vertically) logarithmically spaced grid points; grid 
refinement in pre-assigned rectangular patches, even recursively; 

• parallelization: via MPI (additional OpenMP nearing completion); 

• portability: is running on 

— AMD, IBM and Intel processors 

— AIX, Linux and OSX operating systems. 

The use of high resolution numerics, together with grid refinement, has of 
course the aim of achieving the best possible effective resolution. Given that 
just this is the aim of so many efforts on the observational side, this seems 
to be a worthwhile endeavor also on part of numerical modelling. Regarding 
observations in the solar case, we mention, for example, the Swedish Vacuum 



Telescope or the Sunrise project [11 1 



We have described 2D high resolution simulations of solar granulation in a 
previous paper [3l[. There the main result was the occurrence of strong acous- 
tic pulses generated by the granular downdrafts and moving over considerable 
length, i.e., one, two or even more granule diameters. Furthermore, we observed 
strong vortex patches generated by the integranular downdrafts which remained, 
however, quite stably near the place where they had been generated and in any 
case well below the photosphere. - We should be able to see such phenomena 
in the resolution we achieve in the grid refinement region in the 3D calculation 
which we discuss in the present paper. The question which we discuss concerns 
therefore the nature of 3D solar granulation at high resolution. 
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2. The equations of radiation hydrodynamics (RHD) 

With the ANTARES code, the equations of radiation hydrodynamics are 
solved in one, two, or three dimensions. 

2.1. The hydrodynamical equations 

The equations governing the dynamics of granulation or convection in the 
Sun or in stars (without magnetic activity for the present purpose), are, see e.g. 
(author?) [s^l, the continuity equation 

9tP + V-(pu) = 0, (1) 

the momentum equation 

dt{pu) + V • (puu + pI) = pg + V • r, (2) 

and the total energy equation 

dte + V- (u(e + p)) = p(g • u) + V • (u ■ r) + Q,ad. (3) 

All physical quantities are functions of space and time coordinates, x = [x, y, z)* 
and t, i.e. p =^ p{x,t). The x-direction points in the vertical direction, y- and 
z-direction are the horizontal coordinates. The meaning of the most important 
physical variables used here and in the following is summarized in Table [TJ 



Table 1: Meaning of the main variables 



p 




gas density 


U — 


{u,v,wy 


flow velocity 


pu 




momentum density 


uu 




dyadic product 


p 




gas pressure 


I 




identity matrix 


g = 


(5,0,0)* 


gravitational acceleration 


T 




viscous stress tensor for zero bulk viscosity 

X.y = M {dxjUi + dx.Uj - (V • u)) 






dynamic (molecular) viscosity 


e = 




'total' energy density, i.e. the sum 

of internal and kinetic energy densities 


T 




temperature 


Qiad 


radiative source term 


Xu 




(specific) opacity at frequency v 


K 




radiative conductivity 
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2.2. The radiative transfer (RT) equation 

Any realistic simulation of solar surface (photosphere and upper convection 
zone) flows must include the energy exchange between gas and radiation. This 
process is described by the radiative heating rate Qrad which is an additive term 
in the energy equation. 

In order to determine Qiad the stationary limit of the radiative transfer 
equation 

r-VI^=pxAS^~I,.), (4) 

see (author?) \2§\, is solved for all ray directions r and for all frequencies v 
resulting the specific intensity /i/(r). The stationary limit represented by ^ 
assumes that the temporal variation of the flow is slow and thus its associated 
timescale is long in comparison with that one of photon transfer. Flow veloci- 
ties are required to be non-relativistic (|u| c). Moreover, the radiative energy 
density and its variation as a function of time are assumed to be small in com- 
parison with the total energy density e and its temporal variation. For a general 
discussion of the validity of these equations see (author?) [s^]- - Since local 
thermodynamic equilibrium is assumed to hold and scattering is omitted, the 
source function is the Planck function. Si, ~ By. The (specific) opacity Xv at 
given thermodynamical conditions (p, T) is approximated such that it can be 
interpolated from precomputed tables (see Sect. I4.4.6[) . If the radiation field is 
known, the mean intensity is 

J- = ^Jl-(^)d^ (5) 
and the radiative energy flux is 

= j L{r)vduj (6) 

(see for example (author?) [1^). The radiative heating rate representing the 
difference between absorption and emission can be determined either from 

Qrad = ^T^P j^Xy{Jy - Sv)dv (7) 
or from the equivalent expression 

Q.^d = - j {V ■F,)dv. (8) 

2. 3. The equation of state 

The conservations laws for radiation hydrodynamics ([I])-® are closed by 
the equation of state (EOS) which describes the relation between the thermo- 
dynamic quantities. The EOS used takes into account partial ionization which 
is present in the uppermost few Mm of the Sun. 
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In the simulations presented here reahstic microphysics is included by the 



up-do-date OPAL equation of state [14 1 



For grey radiative transfer Rosseland opacities given by Iglesias and Rogers, 
14 1, extended at low temperature by Alexander and Ferguson, [l[ and (in all 



more recent simulations) by Ferguson et al., [lO|Lare adopted. For non-grey 
radiative transfer parts of the ATLAS 9 package [20[ are used (i.e. opacity dis- 
tribution functions and model atmospheres) to determine bin-averaged opacities 
and source functions. 

For largely arbitrary chemical compositions, i.e. values of X, Y and Z, the 
equation of state is precomputed for a logarithmic mass density-temperature 
grid which provides all required thermodynamical quantities for a given pair 
ip,T). 

Since the system of equations (P)-© is solved in time to obtain mass density 
and total energy density, the internal energy density (total minus kinetic energy 
density) and temperature must be calculated. The latter is interpolated from a 
precomputed grid which provides temperature as a function of logarithmic mass 
density and internal energy density. Alternatively and more precisely, T can be 
obtained by a Newton-Raphson procedure based on the table mentioned in the 
last paragraph. 

Bin-averaged quantities for non-grey radiative transfer are interpolated in a 
precomputed logarithmic mass density vs. pressure grid. 



3. The ANTARES code 

The ANTARES (Advanced Numerical Tool for Astrophysical RESearch) 
code as described here can perform compressible hydrodynamic simulations with 
full radiative transfer for the ID, 2D, and 3D case on a rectangular grid, all 
comprised in one parallelized Fortran90 program. (Spherical coordinates with a 
logically rectangular grid are possible. We concentrate on the straight, equidis- 
tant grid here as it applies to the simulations which we analyze later on in that 
paper.) - Various high-resolution numerical methods (all of the essentially non 
oscillatory, ENO, type, see (author?) Q, (author?) (author?) fl^) are 
implemented. The order of spatial and temporal discretization can be chosen 
arbitrarily to a considerable extent. The viscous terms can be discretized by a 
fourth-order accurate scheme or be replaced by artificial diffusivities (author?) 
(46j . (author?) Q. Optionally, we can also use the Smagorinsky subgrid scale 
model, [i^]. - Furthermore, the ANTARES code allows local grid-refinement in 
rectangular patches, even recursively. Data from the refinement are projected 
up to the coarser level. 

3.1. Structure of the ANTARES code 

Given the physical state q" at time step n a Runge-Kutta step is performed 
to get the physical state q"+^ at time n + 1. Each substep consists of the 
following parts: 
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• Approximation of the spatial derivatives of the hyperbolic part with ENO- 
type schemes fSect. l4?Tt . 

• Discretization of the viscosity: fourth-order discretization of the viscous 
stress tensor or artificial diffusivities (Sect. 

• Determination of the radiative heating rate (Sect. 14. 4p . 

• Evaluation of all other source terms. 

• Update the values for p, pu and e according to the temporal integration 
(Sect.im). 

• Call the equation of state to get all other used thermodynamical quantities 
(Sect. 1231). 

• Apply the boundary conditions for subsequent use (see Sect. 13. 7| ). 

3.2. Numerical grid 

The region of interest D for the presented numerical simulation is a rectangu- 
lar domain. The rectangular computational domain is equipped with a uniform 
numerical grid: gridpoints in x-direction, Ny gridpoints in y-direction and 
Nz gridpoints in z-direction with up to seven ghost cells at each boundary, de- 
pending on the order of the spatial discretization. The mesh sizes are Ax, Ay 
and Az. (In addition, the ANTARES code also offers the option for logarithmic 
grid spacing.) 

The vertical x-range covers x £ [0, Xmax] (a: = is at the top), consisting of 
Nx intervals with length Ax — Xma^/Nx. Similarly, the horizontal ranges cover 
y G [0,ymax] and z e [0, z,nax], consisting of Ny resp. Nz intervals with lengths 
Ay = yma.x/Ny resp. Az = Zma^/N^. 

The grid points (x^, yj, Zk) are the centers of the cells lij.k- For * ~ 1, . . . iV^;, 
j = l,...7Vy and k=l,...Nz 



x^ = (i-^)- Ax, y, = (j - i) • Ay, zj. = (fc - 1) • Az 

and 

Ii,j,k = [x,_i,Xi+i] X [yj_i,yj+i] X [2^.1,2^+1] 

where 



Aa; ^ Ay , Az 

2~ 



x,±i=x^±—, 2/j±i=2/j±— , Zfc±i=zfe± 
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3.3. Grid refinement 

For detailed studies, for example a single granule or the region near a down- 
flowing plume, local grid refinement can be used. First, the region of interest for 
the refined area (dark grey region in Fig. [1} and the (integer) grid refinement 
factor (for each direction) must be specified. The refined area is surrounded 
by ghost cells (light grey region in Fig. [T|) to allow regular interpolation and 
symmetric differentiation near the boundaries. Then, the values of the physical 
quantities at the coarse grid are interpolated to the fine grid and the ghost cells 
(light and dark grey region). This yields the initial state for the simulations. 
Furthermore, a time step for the time evolution at the fine grid Aicfl = ^ 
for an integer N is calculated such that all numerical time step restrictions are 
fulfilled. Typically, Ai is divided by the grid refinement factor in the vertical 
direction unless there are further constraints due to time step restrictions. 

Each step of a simulation with a grid refinement region consists of a Runge- 
Kutta time step on the whole coarse grid. Then, the initial data at the coarse 
grid are interpolated to the ghost cells surrounding the refined area and N 
Runge-Kutta steps are performed on the fine grid. Before each Runge-Kutta 
step or intermediate step, which corresponds to a time ^gr G [^j^ + At], the 
ghost cells are assigned with linear interpolation values from the physical states 
(\{t) and q{t + Ai) from the coarse grid. The values at the grid points of 
the coarse grid which belong to the refined area are overwritten by a volume- 
integrated approximation from the solution (\GR{i + At) at the fine grid after 
N Runge-Kutta steps for the latter (back projection onto the coarser level). 

The back projection step causes a slight violation of conservation. To be 
more precise: in the interior of the refinement region the back projection is 
strictly conservative. However, at the boundaries of the refinement region the 
borders of the refined cells will not necessarily coincide with the borders of 
the coarse cells. There, a non-conservative sort of back projection is presently 
used. No unfavourable effects (such as noticeable mass loss) have been observed. 
Otherwise, the numerics proper is of the conservation type so that mass, mo- 



Figure 1: Domains for a two-dimensional grid refinement. 

mentum and energy are conserved (except for fiow of those quantities through 
the upper or lower boundary). Of course, conservation does also not apply to the 
Qrad term in Eq. [3] in those regions where the radiative transfer equation (and 
not the diffusion approximation) is used, since this term cannot be reasonably 
cast into divergence form in the present context. 
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3.4- Parallelization 

Parallelization is achieved with MPI according to the distributed memory 
concept. Each processor performs a simulation on a rectangular subdomain. 
The number of required processors is specified by the number of subdivisions in 
X-, y-, and z-direction {Px, Py, Pz)- The resulting subdomains are equipped with 
ghost cells in all directions. If an edge area of a subdomain is adjacent to an- 
other subdomain, the ghost cells are filled with the values from the neighboring 
subdomain. Each (intermediate) state is distributed to the adjacent subdomain. 
If an edge area represents a physical boundary then the physical boundary con- 
ditions are applied. Corner values (see Fig. [2]) in any {y x z) -plane (also in the 
ghost regions in x-direction) are interpolated, if required. 

The one-dimensional decomposition procedure is illustrated in Fig. [3) The 
whole region (physical (white) and ghost (grey) cells) is split up. For each 
subdomain additional ghost cells are introduced. These ghost cells are filled 
with the values of the adjacent domain. 

The details of the parallelization of the non-local radiative transfer treatment 
are described in Sect. 14.41 

For the parallelization of a model with a grid refinement zone, the coarse 
and the fine grid are split to the processors according to the distributed memory 
concept. Hence, the computational effort is uniformly distributed among all 
processors. 



z 

— > 



Figure 2: Corners in the (y X 2)-plane. 



3. 5. Performance 

The ANTARES code has been run successfully on up to 256 processors in full 
RHD simulations using MPI parallelization only. Replacing for test purposes 
the radiative heating rate Qiad by the diffusion approximation the algorithm 
scales almost perfectly as long as the subdomains on each processor are large 
enough such that the ghost cells are a small part of the computational domain 
and thus the data transfer between the nodes remains small. 

Since the radiative transfer is non-local, a full RHD simulation cannot scale 
perfectly. Our implementation of the radiative transfer (see Sect. 14. 4p for ex- 
ample yields a theoretical efficiency of 80% (neglecting data transfer and omit- 
ting the transition to the diffusion approximation) for a 3D simulation with 
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distribution of the initial state 
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Figure 3: Domain decomposition in two subdomains for the one-dimensional case. 



Ax = Ay = Az on 64 processors {Px = Py = = 4). In our RHD simulations 
the code is mostly run, however, with P^ equal to 1 or 2, because the diffu- 
sion approximation can be used for the computation of Qrad for those layers of 
the model domain which always remain optically thick. The optimal choice of 
Px, Py, and Pz is thus a compromise that minimizes total computational costs, 
idling due to different physical complexity and hence different equations solved 
within the computational domain (RHD equations vs. diffusion approximation), 
and idling due to data transfer between processors. 

3.6. Initial conditions 

Typically, we use a slightly perturbed static model atmosphere or envelope 
which we equip with a small seed velocity field or pressure perturbations etc. 
when starting a model. For example, for 2D simulations of solar convection 
described in (author?) [sij and furthermore discussed below we start from 
plane parallel models of the solar atmosphere and an envelope due to (author?) 
[6[, with a small velocity field applied. For 3D simulations a converged 2D 
physical state is converted into a 3D one by putting copies of the two-dimensional 
model side-by-side. Again, velocities are perturbed slightly so as to introduce 
a truly 3D field and then the model is evolved in time. Overall, this reduces 
the computational costs for relaxing 3D simulations to a statistically quasi- 
stationary state by up to about 50%. 

3.7. Boundary conditions 

All quantities are assumed to be periodic in both horizontal directions. 

3.7.1. Hydrodynamical boundary conditions 

The solar simulations are currently performed with closed boundary condi- 
tions at the upper and lower boundary of the computational domain. These 
boundary conditions lead to unphysical reflections of waves and shocks. Thus, 
the results are not realistic in the upper photosphere. Presently, our interest 
focusses on lower layers, particularly around and below the superadiabatic layer, 
which essentially are unaffected [19|. 
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For closed boundaries at the top of the computational domain it is assumed 
that the vertical fluxes of mass, horizontal momentum, and internal energy 
vanish. This is achieved by setting the vertical velocity and the vertical gradients 
of mass density and internal energy density to zero there: 

u = 0, dxp = 0, ,dxe = (at the top). (9) 

On the horizontal velocity components stress-free boundary conditions are ap- 
plied: 

dxV = 0, d^w = Q (again at the top). (10) 
At the lower boundary the conditions for momentum read: 

u = 0, dxV = 0, dxW = (bottom). (11) 

If the bottom of the simulation is placed inside a convection zone, as is indeed 
always the case for solar surface convection, the implemented lower boundary 
conditions do not allow fluid motion through this lower boundary. Thus the 
convective and the kinetic energy flux densities arc zero there. In order to 
have the proper energy being transported into the computational domain, we 
introduce a diffusive radiative flux = k*VT. k* is chosen such that the 
introduced diffusive radiative flux approximates the total energy flux at the 
bottom (Fd — » -Ftot)) and such that k* — *■ k, the physical k, three points inside 
the simulation domain (Fd = k*VT converges to F-ad)- This procedure replaces 
the explicit boundary condition on e at the bottom. 

At the top of the domain three ghost cells are used to implement the bound- 
ary conditions. The physical state of the horizontal layer with i = 1 is continued 
vertically to the ghost cells for mass, energy, y- and z-momentum densities. The 
ghost cell value for the a;-momentum density is set to 0. 

At the bottom one ghost cell is used to implement the momentum bound- 
ary conditions, y- and z-momcntum densities are continued vertically, the x- 
momentum density is set to zero. To adjust the radiative flux density the ver- 
tical temperature gradient at each point is constant in time and k is increased 
at three grid points to get the desired radiative flux density. 

No simple rule seems to exist which allows the boundary conditions to be 
implemented in a stable way. So, the exact number of ghost cells actually applied 
to the variables at the top or bottom has emerged from numerous experiments. 

Due to the implemented boundary conditions mass density and energy den- 
sity are conserved. 

High order numerical methods may be in need for more ghost cells than 
we provide. If so, we reduce the order of the method near the boundaries in 
question. Alternatively, we resort to asymmetric stencils. 

3.7.2. Radiative boundary conditions 

To solve the RTE in the computational domain boundary conditions for 
incoming radiation must be specified. 
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All quantities are assumed to be periodic in the horizontal directions. In the 
parallel mode, data transfer between the nodes is performed after each Runge- 
Kutta (sub-)step. Regarding radiative transfer, we need the intensity values at 
the horizontal boundaries of each subdomain. These values are taken from the 
previous (sub-)step. This small time-lag does not noticeably affect the solution. 

At the top of the computational domain it is assumed that there is no in- 
coming radiation, i.e. ^;/(r)ltop = 0- 

At the bottom of the computation domain the diffusion approximation is 
valid, hence /iy(x, r)|j^^j = B^{x). It is valid up to optical depths greater than 
about 10. In all applications we consider the optical depth of the lower boundary 
is orders of magnitude higher than that. 

To save computational time the lower boundary condition is applied at a 
fixed optical depth tda > 10 (calculated with Rosseland mean opacities) and in 
regions with optical depth greater than tda the radiative heating rate Qrad is 
calculated according to the diffusion approximation Qrad = V • (kVT) . 

3.8. Time step restriction 

The maximal allowed time step is determined by the CFL-condition. In a 
single time step the flow is not allowed to transport information by more than 
one mesh width by both motion or sound waves. Thus, we get a time step 
restriction 

A/ <r min,(Aa;,)^ 

^tCFL S L-Courant " T, — HT"; J T (.J-^j 

max(|u|) + max(csound) 

where Ccourant is the Courant number and Cgound is the sound speed. The max 
in the denominator is taken over the whole domain. Ccourant is set to j. 
The viscous stress tensor results in the time step restriction 

A^diffusivc 1^ ^'diffusive T" ' ■ (1^) 

^ max 1/ 

If artificial diffusivities (see Sect. I4.2.ip are used, an additional diffusive time 
step restriction is applied, 

mm^ ( ^^Xi ) 

Aidiffusivo < C'diflusivc ■ (14) 

max J/ 

Cdiffusive is also i. Again the max in the denominator is taken over the whole 
domain. The argument of this max are all shock stabilizing and hyperdiffusive 
viscosity coefficients. If subgrid-scale viscosities were used, time step restrictions 
would occur in a similar manner and thus also accounted for as in equation 1 141 



4. Numerical methods 

The conservation equations (H])-© are treated by the finite volume method. 
The values of the fluxes are interpolated to the cell boundaries and then the 
divergence for each cell is computed. The radiative heating rate is the most 
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expensive part of the whole simulation and is determined by the short char- 
acteristics method for the RTE. All other source terms are simply evaluated 
at the cell centers. Furthermore, the interpolation tool used and the temporal 
discretization are described. 

4.I. Numerical schemes for conservation laws 

The spatial derivatives of the hyperbolic part of the system of conservation 
laws (P), (HI, (131), 

5tq + 9./(q) + 9y.g(q)+a,/i(q) =0, (15) 

where 

u — {u,v,w)* 

q = {p,pu,pv,pw,e)* 
/(q) — {pu, pu'^ + p, puv, puw,u{e + p))* 
g{q) — {pv, pvu, pv'^ + p, pvw,v{e + p))* 
f^-il) — {pw, pwu, pwv, pw'^ + p,'w{e + p))* , 

are calculated for each direction separately. Approximations fi±i j ki 9ij±^k 

and /ij J 1 at the cell boundaries are evaluated in order to get approximations 
for 



9./(q) ^ 






dy9{(\) S 






5./i(q) - 







at each point Xi^^k- 

The values for the numerical flux function at the cell boundaries are cal- 
culated for each component using one of the following high-resolution one- 
dimensional methods: 

• Fourth-order convex non-oscillatory scheme (CNO-4), see 

• Fifth-order weighted essentially non-oscillatory (ENO) scheme (WENO- 
5), see (author?) T^, (author?) Ji^], optionally in combination with 
Marquina flux splitting (author?) [7| 

• Third-order essentially non-oscillatory scheme (ENO-3), cf. (author?) 
0, (author?) 
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4-1.1. Essentially non- oscillatory (ENO) and weighted essentially non- oscillatory 
(WENO) schemes 

Numerical methods for systems of nonlinear conservation laws must capture 
steep gradients (shocks and contact discontinuities) that may spontaneously 
develop and persist in the solution. Classical numerical schemes either produce 
oscillations near steep gradients or smear out these gradients as well as fine 



details, see e.g. [21 1. 



The ENO type methods use adaptive stencils to avoid oscillations in the com- 
puted solution. Typically, when interpolating from cell centers to cell bound- 
aries, around the cell boundary point to which one wants to interpolate the flux 
functions, the stencil consisting of neighbouring cell centers is selected or given 
maximum weight which leads to the smoothest interpoland for the flux functions 
using the stencil points. In addition, upwind considerations are applied. 

For systems of conservation laws this interpolation must be done in the local 
characteristic field since these quantities properly propagate in unique directions 
whereas the changes in density (say) will come about by signals travailing along 
the and the -I- and — characteristics. They arrive, in subsonic flow, from 
both sides and no meaningful upwind philosophy can therefore be devised in 
the original dependent variables. 

Speaking about the x— direction term in Eq. 1151 the eigenvalues and eigen- 
vectors of the Jacobian of flux the function /'(q) at point a;j_|_i ^ (sic!) are 
determined. Then the flux function is transformed for sufficiently many neigh- 
boring points Xi-r:i+s,j,k (again sic!) (r -I- s -|- 1 = fc is the order of the ENO 
approximation) to the local characteristic variables. Then, for each charac- 
teristic variable, the ENO algorithm is performed in order to obtain the cell 
boundary values at x^^i ^ j, for each component. The resulting cell boundary 
values are transformed back. This procedure is repeated similarly for the other 
directions. 

In the ANTARES code a fifth-order accurate weighted ENO scheme and a 
third-order accurate ENO scheme are implemented. For both methods optional 
use of Marquina's flux splitting and entropy fix, see 0] , is possible. 



i.1.2. Convex ENO (CNO) schemes 

Convex ENO schemes as described in (author?) avoid the need to trans- 
form into the characteristic system. They are therefore simpler to implement 
and faster than ordinary ENO methods. They achieve this by splitting the flux 
function / into two components, / = \{f^ + f~). The constitutents are 
designed in such a way that, for /+ all characteristics run from the left to the 
right, and conversely for /~. Thus, for each of them, there is a unique upwind 
direction and no need to go into the characteristic system in order to separate 
the two possible upwind directions. The convex ENO scheme uses, in addition, a 
reference flux, typically a modifled second-order local Lax-Friedrichs flux, which 
reduces to first-order near discontinuities, and chooses the convex combination 
of the interpolation values of the different candidates stencils which is nearest 
to the reference flux. By the choice of a reference flux depending smoothly on 
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the parameters, albeit of low order, Lax-Friedrichs in our case, and by this pro- 
cedure smoothness of the resulting flux is achieved in addition to its high order 
due to the use of high order polynomials for interpolation. 

More specifically, we consider now the evolution equation for one of our basic 
dependent variables and describe the procedure obtaining cell boundary fluxes 
in direction x (index i). For the various approximations (due to various choice 
of the stencils) at the cell boundaries for a component of the two flux functions 
for one component of our conservative system of equations 



and 

where, in principle. 



fti,j<i) = km+'^.+i,,k<i) (16) 



f-+i,,,i<i) = lim-^.+uk<i) (17) 



max \f'{q)\ (18) 

min{qij^k ■qi+l,i,k)<q<mayi{qi^j_k,qi+l,j,k) 



given by the ENO algorithm in x-direction, the CNO algorithm in the ANTARES 
code chooses for each sign ± the flux which is nearest to the Lax-Friedrichs ref- 
erence flux, if all candidate fluxes are all larger or all smaller than the reference 
flux. Otherwise, the reference flux is taken. 

The choice of a^^i ^ ^. as in equation [T51 guarantees that, for each of the 
auxiliary fluxes the characteristics point all into one direction. In practices, 
a can be and is taken somewhat larger than the value just quoted resulting in 
greater stability at the expense of more smearing. 

4.I.3. Test problem 

To show the main difference between these three schemes the one-dimensional 
version of the system ([T5|) with the initial conditions 



\ = \ I forx<^ (19) 
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for X > ^ (20) 

at 60 spatial grid points is evolved 100 second-order Runge-Kutta steps with 
At =0.7s. The resulting state is very similar for all this three high-resolution 
methods (see Fig. S]). 

A closer look at regions with steep gradients (see Fig. [5]) shows some differ- 
ences between the three methods. The fifth-order WENO scheme produces the 
steepest gradients. The CNO scheme is the most smoothing one. The third- 
order ENO solution lies between the two others in regions of steep gradients. 
The difference in the solution between the two different ENO-type schemes is 
due to the order of accuracy. 
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4-2. Viscous hydrodynamical fluxes 

In addition to the numerical diffusion of the ENO/CNO-schemes for the 
conservation laws the physical viscosity can be included by a fourth-order dis- 
cretization of the viscous stress tensor or it can be replaced by artificial diffu- 
sivities. 

4-2.1. Discretization of viscous fluxes 

The derivatives in the viscous stress tensor t are calculated by the fourth- 
order scheme 

du\ 1 

Then, the values for r and (u-t) at the cell centers in the momentum and energy 
equations are interpolated to the cell boundaries by the following scheme: 



and the equivalent is done for f^_i to get the contribution by V-r and V- (u-r) 
through 



Ax 

For simulations with this implementation of the viscous stress tensor we 
can choose a Prandtl number. Note that the Prandtl number typical for the 
solar atmosphere, see Pr = varies from 1.4 • 10^'* at the bottom to 
1.0 • 10~^° at the top of the computational domain. For our actual choice see 
the discussion at the end of subsection l5.4l - Wc note that the radiative Prandtl 
number defined via /irad = pt^rad = 3Tk/(4c^), cf. 16], is generally larger in 
the layers near the photosphere than the molecular viscosity; the latter has e.g. 
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Figure 5: Comparison of the numerical schemes for a region with steep gradients. 



been calculated in For our simulations of solar granulation Pr^ad ranges 
from a few times 10~^° to about 2 • 10~^. As Pri-ad can easily be calculated from 
equation of state related quantities and opacities which are already required to 
calculate the radiative energy transfer, only the radiative viscosity is used in 
the simulations presented here. Hence, our Pr should more accurately be called 
Prrad, since conductive contributions to k are negligible for the solar surface 
layers. It is important to not confuse this definition of Pr with a purely molecular 
Prandtl number that compares molecular viscosity with heat conductivity: in 
most regions inside stars, typically as long as the equation of state describes 
a perfect, non-degenerate gas, the heat conductivity is orders of magnitudes 
smaller than the radiative one. Indeed, in stars the Prandtl number defined 
through Prcond = TT^^ much larger than Prrad and much closer to that one 
of air. Physically, Prcond is irrelevant, because either radiation or convection 
or both outperform conduction in heat transport inside the Sun by orders of 
magnitudes. 

Finally, we note that in those cases where we only use a physical viscosity 
term instead of one that is enhanced by an artificial diffusivity model, the nu- 
merical viscosity of the method for the hyperbolic part must be able to stabilize 
the simulation. 

4-2.2. Artificial dijfusivities 

Artificial diffusivities as described in and remove short-wavelength 
noise without damping longer wavelengths and diffuse strong discontinuities in 
order to stabilize the numerical code. 

The viscous stress tensor 




(21) 
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is replaced by artificial equivalents of the form 

Zki = ^pMu)d^,Uk + i'i{u)d^^ui) (22) 

where u = (iti, U2, u^)* and k,l — 1, 2, 3. 

The coefficients I'k for direction k consist of two parts, a shock resolving i/f^^ 
and a hyperdiffusive i^]^^^- 

^fc(u) = ^f'^(u)+z.,^''PK). (23) 

The shock resolving part in direction k is defined by 

shk. ,x / C^shkA4|V-u| V-u<0 
W-<^ Q V • u > 

and acts therefore in regions undergoing compression only, i.e. where V • u < 0. 
The hyperdiffusive part in direction k is defined by 

hyp.n A max3 \Alf\ 

I'k if) = ChypAx,ctot (25) 
max3 I A J. /I 

Thereby we have 

Ctot — 1^1 ^" Cgound; 

A\f denotes the first difference of a scalar quantity / (uk in the actual ap- 
plication) in direction fc; maxa indicates that the maximum is taken over three 
adjacent intervals. A^/ designates a similar third difference. For all simulations 
presented here Cghk = 1 and Chyp — 0.05. These values make reasonable contact 
with those used in [HI, namely Cghk = 2 and Chyp = 0.05, and in fHo*!, Cghk — 1 
and Chyp = 0.03. 

In addition, we have implemented the Smagorinsky subgrid model, see (au- 
thor?) [31, for optional use. 

4-3. Interpolation tools 

We need two interpolation schemes: one for use in the equation of state tables 
(in order to get smooth derivatives) and in connection with grid refinement, 
and an other one, being applied in the radiative transfer solver (in order to get 
positive opacities etc. even in the presence of huge gradients), which we describe 
in turn. 

For EOS and grid refinement interpolation we use local splines [3^ . Starting 
with four points, two polynomials through both of the three adjacent points are 
used to get two interpolation values at point P. These two results are combined 
by a weighted sum. The weight for polynomial 1 with stencil 1 is ^^'^^^ and for 

polynomial 2 with stencil 2 is (see Fig. [B]). In the two-dimensional case 

4x4 grid points are used to get the interpolation values. First, four one-dimen- 
sional interpolations in direction 1 and then one one-dimensional interpolation 
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in direction 2 are performed to get the value at point P. The successive inter- 
polation procedure is illustrated in Fig. [T] For three-dimensional interpolations 
4x4x4 grid points are used and the procedure described above is performed 
starting with 16 interpolations in direction 1, followed by 4 interpolations in di- 
rection 2 and then one interpolation in direction 3. This method yields smooth 
(continuous) first derivatives everywher (cf. [s^ for details on an efficient im- 
plementation also applicable to vector fields). 

In the radiative transfer solver a method for monotonic interpolation in one 
dimension proposed by StefFen is applied. The method gives exact results 
if the data points correspond to a second-order polynomial. The application 
of this methods assures that the values for density, opacity, intensity, and the 
radiative source function always remain positive. 

Stencil 2 




Stencil 1 

Figure 6: Stencils for one-dimensional interpolation. 



P 

O 



Figure 7: Two-dimensional interpolation procedure. 



4- 4- Numerical radiative transfer 

Instead of solving the three-dimensional RTF ^ the one-dimensional RTF 



(26) 



is solved along several ray directions, following the short characteristic method 
proposed by Mihalas [2^ and Kunasz and Auer [3]. dr^ = XuPdx is the optical 
depth of a path element at frequency v along a ray direction. 

In the following the subscript v which indicates the frequency dependence 
of /, J, S and r is dropped. Thus, the described numerical schemes treats grey 
radiative transfer. In Sect. I4.4TB1 frequency dependent RT is discussed. 
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4- 4-1- Short characteristic method 

For incoming radiation the values for I{t[U)) along ray direction r are 
known. The optical depth t{U) there is set to zero. To determine I{t{P)) 
along this ray direction r (see Fig. [S]) one interpolates the values of the required 
physical quantities to the point U using 16 points in 3D and 4 points in 2D, 
respectively. Then, one evaluates the equation 

/(r(P)) = /(T(C/))e^(^)-^(^) + 

/ 5(r')e" -"(^^dr' (27) 

Jt{U) 

numerically to get I{t[P)). This procedure is repeated recursively, since after 
step 1 the intensities for the first layer are determined (see Fig. [5] and ^ . 
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Figure 8: Short characteristic method, step 1 for ray directions entering at the top. 
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Figure 9; Short characteristic method, step 2 for ray directions entering at the top. 

4-4-^- Numerical integration 

To perform the numerical integration on the right hand side of (j27p one has to 
determine t{P) and t{P + UP), where t{U) is set to zero. The approximations 
are calculated by twice evaluating J pxdx using the trapezoidal rule. 

The integration on the right hand side of Eq. [27]is performed by a quadrature 
rule proposed in [s^. 
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4.4-3. Angular integration 

Table 2: Support points Xi and the weights for the Gaussian quadrature rule for A'^ = 10 
(see [H) 



n=10 i 




at 


1 


-0.974 


0.067 


2 


-0.865 


0.150 


3 


-0.680 


0.219 


4 


-0.433 


0.270 


5 


-0.149 


0.296 


6 


0.149 


0.296 


7 


0.433 


0.270 


8 


0.680 


0.219 


9 


0.865 


0.150 


10 


0.974 


0.067 



For every grid point Xij_k the RTE is solved along iVj-ays ray directions. To 
get the mean intensity J an angular integration is performed. 

In the two-dimensional case the RTE is solved along 20 ray directions = 
(■Til, ri2, = 0). The ra are set equal to the support points Xi of the Gaussian 



quadrature rule for = 10 according to [2J], the ri2 are determined such that 
'"'ii ~^''^i2 ~ ^ e^or every rn there are two solutions for rn, one with positive, one 
with negative sign). For the angular integration the corresponding weights 
are multiplied by tt since 



10 

E 

4=1 



= 4tt (28) 



is required. The value of An instead of the possibly expected 27r stems from the 
fact what we consider to be the proper transport of radiation in 2D keeping the 
3D case in mind as the ultimate goal. Namely, at a fixed location we do not only 
consider one ray of light (corresponding to a, say, positive horizontal direction 
and a certain zenith distance) within the computational plane but rather the 
half cone generated by such rays, which are allowed to move out of the plane 
of computation, pointing into positive horizontal direction and the same zenith 
distance). - The values for Xi and are shown in Table [21 

In addition to this standard integration, ANTARES also offers the possibility 
to use 12 rays in 2D and to use the weights — 2/(N_rays/2). 

In the three-dimensional case the RTE is solved along 24 ray directions. 
The ray directions are chosen according to the angular quadrature formulae 
of type A of The directions in each octant are arranged in a triangular 
pattern and the quadrature is invariant under rotations over multiples of 7r/2 
around any coordinate axis. A summary of the construction procedure is given 
m (author?) ^. For the three-dimensional simulations the A4 quadrature set 
with three directions per octant is used. 
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For this choice each ray has got the weight uji = The ray directions in 
the first octant are: 



= (v/779,1/3,1/3) 
= (1/3, ^779,1/3) 
= (1/3, 1/3, V779) 



The ray directions in the other octants are obtained applying three-dimensional 
rotations into the other octants. 

With the quadrature weights tOi and the index i running over the set of 
directions, the mean intensity is 



4- 4-4- Radiative heating rate 

With J now calculated at each grid point the radiative heating rate can be 
obtained from 



4- 4- 5. Implementation 

Radiative transfer is a non-local phenomenon. Simulations with Np pro- 
cessors and MPI communication are carried out by the domain decomposition 
approach. Since the intensity is calculated with the short characteristic method 
along A^iays rays which start at a specific boundary of the computational do- 
main, each processor has to wait until the information at the specific adjacent 
boundary is available. 

For simulations with numerous processors one has to ensure that as many 
processors as possible are busy. Our short characteristic algorithm classifies the 
A'lays rays, depending on the direction in which they pass the computational 
domain. 

Considering the short characteristic algorithm, the rays are numbered such 
that the iVj-ays.x/S rays entering at the top of the cell are followed by the Arrays, x/2 
rays which enter at the bottom, followed by the Arrays, y/2 rays entering at the 
lower y-boundary and by the iViays.y/S rays entering at the upper y-boundary, 
followed by the N^a.ys,^/'^ rays entering at the lower z-boundary and by the 
^rays,z/2 rays entering at the upper z-boundary. 

For an equidistant three-dimensional grid the four ray directions entering at 
the top of the cell are called r^, . . . r"*. The four ray directions entering at the 
bottom of the cell are called . . . r^. Of course, for the remaining four faces of 
the computational domain there exist similar directions r^, . . . r^^. 

Figure \W\ illustrates the implementation for the ray directions entering at 
the top and the bottom of the cell. The vertical axis is devided in P^, = 6 parts 
which represent the 6 subdomains in x-direction. 




(29) 



Qrad = 4:T^PX{J - S). 



(30) 
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The processors representing the upper half of the computational domain first 
calculate the intensities for the rays entering at the top of the cell (r^, . . .r'*) 
in the following way: In the first step the processors representing the first layer 
from the top calculate the intensity for r^. In the second step these processors 
calculate the intensity for and the processors representing the second layer 
from the top calculate the intensity for r^. This procedure is continued until 
all processors representing the upper half of the computational domain have 
calculated the intensities for all rays entering at the top of the cell. 

The processors representing the lower half of the computational domain first 
calculate the intensities for the rays entering at the bottom of the cell (r^, . . . r®). 
The procedure is similar to the one described above, only the start is at the 
processors representing the first layer from the bottom and then moving forward 
to the top. 

After that, the processors representing the upper half of the computational 
domain have to continue from the bottom to the top with the rays entering at 
the bottom of the cell (r^, . . . r^) and the processors representing the lower half 
of the computational domain have to continue from the top to the bottom with 
the rays entering at the top of the cell (r^, . . . r"*) in the way described above. 

The efficiency of this algorithm depends on the relation between P^. and 
-^rays,x- The morc subdomains in x-direction and the less ray directions entering 
at the top or at the bottom, the less efficient is the implementation. 

The algorithm is performed componentwise, first the rays entering at the 
a;-borders, then the rays entering at the y-borders, and finally the rays entering 
at the z-borders are considered. 



g step of the algorittim 

'm I ► 

E 
o 



1 


2 


3 


4 










5 


6 


7 


8 




1 


2 


3 


4 






5 


6 


7 


8 








1 


2 


3 


4 


5 


6 


7 


8 










5 


6 


7 


8 


1 


2 


3 


4 








5 


6 


7 


8 






1 


2 


3 


4 




5 


6 


7 


8 










1 


2 


3 


4 



Figure 10: Illustration of the implementation of the radiative transfer solver. Rays with 
number 1, 2, 3, 4 enter at the top of the cell, rays with number 5, 6, 7, 8 enter at the bottom. 



4.4-S. Frequency dependent radiative transfer 

The frequency dependent radiative transfer, which is important in the solar 



photosphere, is implemented by the opacity binning method, see [35[ and |25j . 

Frequencies v which reach optical depth Ti, = 1 at similar geometric depths 
in a one-dimcnsional atmosphere are grouped into the same bin. It is assumed 
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that frequencies which reach optical depth unity at similar geometrical depth 
have got a similar behavior in the radiation field. 

For each frequency set fti a mean opacity Xi ^^'^ ^ mean source function 
Si are determined for which the radiative transfer equation is solved, resulting 
in Ii{r). Then the mean intensity Ji for each bin is calculated. Finally, the 
radiative heating rate is given by the weighted sum over the bins, 

Qrad = 47rp ^ U!iXi{Ji " Si). (31) 

1=1 

In the ANTARES code 4 or 12 bins are used. Grey simulations are carried out 
with the Rosseland mean opacity and the frequency-integrated Planck-function. 
The procedure used for calculating the binned opacities is based on a pack- 



age for computing opacity distribution functions by [37|. Though in principle 
we could also use this package to compute new opacities for arbitrary abun- 
dances, in the present work we start from the precomputed tables of Kurucz, 
(20| . as we first want to study the possible improvements gained from a more 
refined numerical treatment of the hydrodynamics. The software package of 
Piskunov & Kupka, |37|. also includes a code for adding continuum opacities to 
the tables of Kurucz, [20|, which only contain line opacities. Summation over 
the resulting opacity distribution yields Rosseland and Planck opacities. This 
capability was required since they computed opacity distribution functions for 
chemically peculiar stars for which the precomputed Rosseland tables [2^ could 
not be used. 

To compute binned opacities we take a one-dimensional model atmosphere 
as in [20| — or in principle also horizontal and temporal averages from an ear- 
lier numerical simulation with the same basic parameters — to determine the 
Rosseland optical depth tross at which each bin of the opacity distribution func- 
tion reaches Ty(TR,oss) = 1- Following this analysis each bin is assigned to its 
corresponding optical depth bin. On this basis we calculate partial sums of 
frequencies to which these bins correspond as well as of the Planck function 
and its temperature derivative, in addition to Rosseland and Planck means. 
By performing a summation of these quantities over all frequencies of interest 
(from about 9 nm to 160 /im) we obtain the frequency weight, Rosseland and 
Planck means as well the partially summed Planck function and its derivative 
with temperature for each optical depth bin. Rosseland and Planck opacities 
for the entire frequency range are computed on the fly and compared with a 
computation based on integrating the optical depth bins in the same manner 
(harmonic mean weighted with the temperature derivative of the Planck func- 
tion for the former and an arithmetic mean weighted with the Planck function 
but excluding contributions from scattering for the latter). For the tables used 
in our simulations these averages agree reasonably well. For the case of Rosse- 
land opacities we have also successfully compared our results to the original 
tables of Kurucz, [2^, as was done already in ^ for the case of A and B type 
stars. Overall, our procedure is very similar to that one discussed in detail in 
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For the opacities eventually assigned to each optical depth bin i we go 
one step beyond the procedure just described and follow [25| by performing 
a weighted average of Planck and Rosseland opacity as a function of depth. 
Thus, 

% = e"""^"--' Xpianck,* + (1 - C"^"°--") XRoss.i (32) 

where trqss is estimated from a simplified hydrostatic equilibrium relation: 
TRoss.i ~ Xross iP/d- This puts the transition region between Rosseland opacities 
(high optical depth) and Planck opacities around tross « 1/e. This setting has 
been used in all the calculations shown here, but our binning code also allows 
changing the location of this transition. In principle, this parameter can be 
optimized together with the number of bins and their inner boundaries. 

In our simulations for grey radiative transfer we use the Rosseland mean 
computed from these tables. For the case of simulations in two spatial dimen- 
sions and non-grey radiative transfer we mostly use 12 bins delimited by the 
set of points logTRo^s = -hO.25, 0.0, -0.25, -0.5,-1.0,-1.5,-2.0,-2.5,-3.0, 
—4.0, —5.0, while for simulations in three spatial dimensions we mostly use four 
bins delimited by logrRoss — -0.5,-1.5,-2.5, as in [2^ for the case of cool 
stars. 

4-5. Temporal integration 

To advance the numerical solution in time either a second or a third order 
accurate Runge-Kutta scheme, (author?) [i^l, is used. These time integra- 
tion schemes are specifically designed for use in conjunction with ENO spatial 
discretisations with the aim of not introducing unwanted additional spatial vari- 
ation. 

5. Comparison of the numerical methods 

In order to get insight on the quality of various numerical methods, we cal- 
culate two-dimensional models of solar granulation, physically identical, with 
different numerical schemes. Due to the dimensional splitting approach the nu- 
merics in two and three dimensions is as similar at possible at all. Hence it 
seems reasonable that relative advantages and disadvantages of various numer- 
ical schemes carry over from 2D, where such experiments are cheaper, to 3D. - 
In addition, we compare two 3D models. 

5.1. Simulations 

Here we discuss the influence of numerical methods for the hyperbolic part 
(WENO-5, ENO-3 and CNO-4) including artificial diffusivities (AD) with WENO- 
5 and ENO-3 methods without AD to investigate the influence of the order of 
the methods used for the hyperbolic part of the entire problem. We also consider 
the influence of AD per se (WENO-5 with AD and WENO-5 with a minimal 
Prandtl number) . The Prandtl number Pr is an important similarity parameter 
for fluid flow: 

Pr=^, (33) 
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jj, is the dynamic viscosity, k is the thermal (radiative) conductivity and Cp is the 
specific heat at constant pressure. The right hand side of ((33)) can be rewritten 
as t^/'' ^ . Hence, the Prandtl number is the ratio of the kinematic viscosity to 
the thermal (in our case actually radiative) diffusivity. 

Through the comparison of the entries t^,; of the viscous tensor, 

Pr^ (^d^^Uk + dx.ui - ^4;(V • u)^ (34) 

and the corresponding component of the artificial viscosity, 

]j^p{vi{uk)dxiUk + Vk{ui)dxt,ui) (35) 

one can estimate a formal "Prandtl number" Pr which in fact is rather a Peclet 
number defined for the scale of the artificial diffusivities: Prad Pcad = 
Vad/ {n/ {cpp)). Figure [TT] shows this estimate for the initial configuration from 
which we begin our further comparisons of the two-dimensional numerical sim- 
ulations. At the top this formal Prandtl number is a few times 10~^, increases 
rapidly and attains a value of 10^ at the bottom. This demonstrates that near 
the top a negligible amount of heat is transported at the grid scale and con- 
firms that for the surface layers the physically important contributions to energy 
transfer are resolved on the computational grid. However, since our simulations 
have to neglect several decades of scales down to those of viscous dissipation, 
once we consider layers deep inside the convection zone, even at the length scale 
of the computational grid convection is found to be much more efficient in trans- 
porting heat than radiation and we expect a huge increase of the formal Prandtl 
number as is indeed found in Fig. [TlJ From Eq. ([34]) we see that for constant Pr 
the contribution by the viscous stress tensor r^,; declines with increasing depth, 
since the fraction n/cp declines rapidly and the fields of the velocity gradients 
do not vary significantly inside the computational domain. On the other hand, 
the entries in the artificial diffusion tensor, Eq. (|35p . increase, since the value in 
the brackets, which shows the same behaviour inside the computational domain, 
is multiplied by the mass density. Since the Prandtl number in the solar photo- 
sphere is about 10"^*^, see subsection 14 .2. 11 the artificial diffusivities serve only 
stabilize the code and do by no means describe the physical viscous behavior. 
Therefore, a comparison of simulations with and without artificial diffusivities 
is essential. As the numerical viscosity of the WENO-5 and ENO-3 schemes 
themselves is difficult to compute and it is, moreover, not clear whether it can 
at all be reasonably described by viscosity terms of the usual form, we resort to 
a comparison of simulation results and their dependency from the order of the 
numerical method used and the presence or absence of AD in the simulation. 

5.2. Simulation setup 

The vertical and horizontal extents of the computational domain are 2763 km 
X 11172 km. This domain is provided with 182 x 485 cells which yields a 
mesh size of 15.2 km in the vertical and 23.0 km in the horizontal direction. 
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Figure 11: Estimate for logjQ of the formal Prandtl number defined from comparing the 
artificial diffusivities with the viscous stress tensor. 



This resolution is common for simulations of the surface layers of the Sun. 
Approximately 80 minutes of solar granulation are simulated. 

5.3. Results 

The temporal means of the horizontally averaged temperature show small 
differences which originate from the advection schemes used for each simulation 
(Fig. [12] displays the region containing the superadiabatic layer). Among the 
simulations with artificial diffusivities, the ENO-based schemes produce steeper 
gradients, which agrees with the observations made for the test problem in 
Sect. 14.131 This effect is obvious when comparing the solutions for EN0-3-AD 
and CN0-4-AD. The simulation without artificial diffusivities yields a different 
temperature profile: a steeper gradient in the upper part (at the solar surface) 
and a flatter gradient in the lower part (underneath it). 

The temporal means of the horizontally averaged velocity components and 
their spatial derivatives (see Fig. [Olandll4lfor two examples) are similar for all 
three methods using artificial diffusivities. The shapes of these curves are all 
similar, though their magnitudes are slightly different. Although the simulation 
without artificial diffusivities also yields similar averaged velocity profiles, the 
numerical values obtained from this method are clearly higher. It both yields 
greater shear stresses and velocity divergences (Fig. [T5)) as well as higher vertical 
velocities (Fig. [Tl)) . 

At a depth of 600 km the granulation simulation is comparable with solar 
observations for determining the mean size of the granules. The mean distance 
d between granular centers, the cell size, was measured to be 1.94" in (4ol |. 
while in a value of d = 1.76" is reported, which corresponds to 1400 km 
and 1270 km, respectively. Although the horizontally averaged quantities are 
similar for the simulations without artificial diffusivities, there is a difference 
in the mean number of downfiows iVoF and the resulting mean cell size d at 
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Figure 12: Vertical temperature profile. Note the shift of up to about 30 km, or roughly two 
grid cells, within the region where the temperature gradient is steepest. 



a depth of 600 km (see Table [3]) . For the chosen resolution none of the three 
methods using artificial diffusivities yields an appropriate number of downflows. 
Only for the simulation with a plain fourth-order discretization of the viscous 
stress tensor (in addition to its intrinsic, small numerical diffusivity) one obtains 
an acceptable mean number of downflows and a corresponding granule size. In 
evaluating this result one should keep the 2D nature of the simulations in mind, 
which may not, in such detail, be comparable to observations. However, our 3D 
simulations result also in smaller granules for even finer grid spacing. Hence, 
these results may also be taken in favor of WENO-5 without AD. 



Table 3: Number of downflows A'^df B.nd granule size for the 2D simulations using different 
advection schemes. 





Nop 


d [km] 


WENO-5 


8.1 


1380 


WENO-5- AD 


7.1 


1580 


EN0-3-AD 


6.2 


1810 


CN0-4-AD 


6.8 


1650 



In addition to one-dimensional averages and the size of the granules we can 
also consider the pressure distribution (logarithm of pressure with logarithm 
of its horizontal average subtracted, see Fig. [T5)) . The simulations with this 
typical resolution for solar granulation simulations of 15.2 km vertically and 
23 km horizontally can be compared with the high-resolution model described 
in ^31j . The acoustic pulses which are ubiquitous in the high-resolution model 
are present in the simulation without artificial diffusivities. In the simulations 
with artificial diffusivities these acoustic pulses are not that sharp and numerous 
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Figure 13: The quantity j (dxjUi) j is a measure for both the turbulent behaviour 

and the presence of shock fronts in the system, since it is sensitive to steep gradients caused 
by shear and by strong compression or expansion flows. 



and not found so deeply down. Furthermore, the simulations without artificial 
diffusivitics show more small scale features (e.g. at the upper right corner) and 
more whirls which we have also observed in our higher resolution runs. 

5. 4-. Discussion 

For simulations with higher resolution and artificial diffusivitics the values 
for the horizonal averages increase and the mean size of the granules decreases 
such that the differences caused by various implementations of the viscosity (or, 
from a different point of view, of the modeling of the unresolved scales) vanish. 
A 2D high-resolution model of solar granulation with artificial diffusivitics in [3l[ 
shows agreement on an quantitative level: the size of the granules is appropriate. 

For 3D simulations we mostly use the fifth-order WENO scheme with Mar- 
quina's flux splitting and entropy fix. This method is stable for both implemen- 
tations of the viscosity (numerical viscosity plus Navier-Stokes type viscosity 
only, or with artificial diffusivitics included) and the mean size of the granules 
is nearest to the observed values for the 2D simulations already at a typical 
resolution for solar granulation calculations of about 20 km. Moreover, if we 
consider a horizontal profile, where the density shows one smooth peak and at 
this peak the vertical velocity changes sign, we find that the temporal develop- 
ment of the density calculated with Marquina's flux splitting and the entropy 
fix remains smooth, whereas the solution calculated with the standard ENO- 
or WENO-method generates oscillations in the smooth peak of the density. 
These observations also correspond to the theoretical considerations described 
in [§]. Thus, all subsequent simulations use WENO-5 and include Marquina's 
flux splitting and entropy fix. 
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Figure 14: The quantity = u^, the squared vertical velocity, is found to be slightly higher 
in the simulation without artificial diffusivity. Note that there is no simple connection to 
kinetic energy flux as the latter depends on flow assymetry. 



If we repeat some of our numerical experiments with these settings in 3D, we 
again find that a simulation without artificial diffusivities using the same mesh 
size yields a markedly higher resolution than its counterpart using artificial 
diffusivities. This is evident from the horizontal slices at a depth of 600 km 
for 3D simulations shown in Fig. 1161 where more small scale features become 
visible. For those two simulations the computational domain comprises 2763 km 
X 11172 km X 11172 km and is provided with 182 x 285 x 285 ceUs which yields 
a cell size of 15.2 km x 39.2 km x 39.2 km. Given the veracity of the methods 
without artificial diffusivity as discussed in the 2D case above we feel justified 
to trust the results in 3D to a high degree as well. 

If we apply WENO-5 with no artificial diffusivities in 3D we have to add a 
Navier-Stokes-type diffusion term in the momentum equation with a coefficient 
larger than the physical one in order to keep the code stable. Considering an ever 
better resolved series of calculations, we can drastically reduce the amount of 
viscosity we need to add as we increase resolution. In fact, in the refinement zone 
of the 3D model discussed in the next chapter we can do by adding a viscous term 
(of the form of usual hydrodynamic viscosity) to the momentum equation only, 
the coefficient being that small that, for the corresponding Prandtl number, 
we have Pr < 10~^ everywhere. Keeping in mind the adverse effects of added 
viscosity as discussed earlier in the 2D case this implies that gain in effective 
resolution is larger when moving to a highly resolved simulations than can be 
expected on basis of the grid sizes alone and that, at high resolution, WENO-5 
can do with, by any numerical standards, only a very small amount of diffusivity 
added. From a somewhat different standpoint, it is remarkable in itself that 
WENO-5 really feels a viscosity which is so small in terms of Prandtl number. 
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6. A three-dimensional simulation with local grid refinement 



We just have mentioned the model with the basic cube of size 2763 km x 
11172 km X 11172 km (equipped with 182 x 285 x 285 cells, cell size 15.2 km 
X 39.2 km x 39.2 km). Into that domain we have plugged a refinement zone of 
size 1966 km x 4165 km x 3851 km, provided with 259 x 425 x 393 cells. That 
leads to a cell size of 7.1 km x 9.8 km x 9.8 km. The high-resolution domain 
starts quite high in the atmosphere. We now discuss the results borne out from 
this high resolution simulation. Basically, our simulation follows the decay of a 
large exploding granule occupying the largest part of the refinement region. 

Quite soon, due to the higher resolution, new features become visible which 
were at best hinted in the lower resolution run. In fact, a considerable number of 
vortex tubes commenced evolving. They are generated by the downflow lanes 
at the boundaries of the convection cells at a depth where the downfiowing 
material gets buoyant and is defiected sidewards and ultimately upwards. 

A number of these tornadoes manages to move upwards into the photo- 
sphere. They can readily be seen in Fig. 1171 The figure contains an isosurface 
plot of temperature, T — 6000 K. The tornadoes are visualized by means of a 
volume rendering of the norm of vorticity. It is obvious that they crowd pre- 
dominantly in the vicinity of the downflows, characterized by depressions in 
the T— isosurface, whereas the broad upflows are void of them. The tornadoes 
can also readily be recognized in other variables: V(p — p) where p denotes the 
horizontally averaged pressure, V(/5 — p), etc. 

The effects are large. Pressure inside the tornadoes may amount to only ^ ^ 
of the ambient pressure, density may also be quite low and, in addition, tem- 
perature is lower than in the surrounding medium. The tornadoes are stabilized 
against collapse resulting from the higher ambient pressure by rapid rotation. 

The vortex tubes reach the surface layers by ascending predominantly at 
the granule boundaries where the vertical velocity is and there exists strong 
shear. Typical diameters are somewhat above 100 km and smaller. Note that 
a tornado of that size has only about 10 grid points to represent it. Given the 
strong velocity gradients it is remarkable that the tornadoes survive. It would 
be an interesting question to study the stability with either analytical models 
or numerically applying even higher resolution. 

These tornadoes can be considered to be the 3D analogs of the vorticity 
patches described, for the 2D case, by ^ij. However, in moving up to the 
surface they behave quite differently than their 2D counterparts which stay 
deep below the photosphere. In fact, in 3D the tornadoes can reach up to the 
top of the computational domain, quite high in the photosphere, as witnessed 
by the strong vortex tube displayed in Fig [T8l 

Obviously, the tornadoes are candidates to explain the extra turbulence near 
downdrafts which has been diagnosed by high resolution spectroscopy of solar 
granulation; see e.g. [33|,[l3|- As it is of interest to investigate those phenomena 
in ordinary granules, we have presently such work ongoing. 

In the paper on high resolution 2D models of solar granulation, 3l|, we 
have observed strong acoustic pulses. Actually, they are also easily seen in 
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Fig. [13] of the present paper which refer to the 2D case, which come from a 
model of much lower resolution than in the paper just mentioned. As much as 
resolution is being concerned their 3D counterparts should clearly show up in the 
present refined 3D resolution. There are, however, no such isolated and strong 
acoustic pulses. Rather, wavetrains can be observed (best in movies where one 
can actually follow them). It will require a more subtle analysis to figure out 
whether their energetics is comparable to that of the pulses in 2D and whether 
they possibly play a role for the heating of the chromosphere. 

7. Summary and conclusions 

We have described the ANTARES code for stellar radiation hydrodynamics 
whose design goals have been flexibility regarding numerical schemes and dimen- 
sions of the computation as well as microphysics (idealized or realistic). That 
holds also true for grid structure (straight or polar coordinates) and equations 
used (e.g. hydrodynamic equations in straight or polar coordinates). 

Then, we have compared various numerical schemes and demonstrated the 
benefits of using high-order ENO schemes. In the high-resolution runs, WENO-5 
with Marquina flux splitting stably with but a tiny amount of viscosity added in 
the momentum equation and without a need to add any unphysical diffusitivity 
in the energy equation. This results in high information content per grid point. 
Grid refinement, possibly in a hierarchical fashion, allows to closer zoom in on 
interesting features. 

In this way, we have investigated solar granulation, in particular a decaying 
exploding granule in the refinement region. Numerous rapidly spinning vortex 
tubes (diameter about 100 km or less) generated by the granule border down- 
fiows in the region where they become buoyant again manage to ascend to the 
photosphere near the granular downflows. Many of them are arclike, some others 
reach straight up to the top of the domain, i.e. high into the solar atmosphere 
and possibly therefore the chromosphere. 

Applications of the ANTARES code to other topics in stellar physics are 
underway and will be reported in due time. 
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(a) WENO-5 




(c) CNO-4-AD 



Figure 15: Isolines of perturbations in the pressure distribution (logarithm of pressure with 
logarithm of its horizontal average subtracted) after 25 minutes of evolution from a thermally 
relaxed model. The isolines crowd near shocks and acoustic pulses. 
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(b) WENO-5-AD 



Figure 16: Horizontal slice of the temperature distribution at a depth of 600 km seven minutes 
after the simulation began with the same initial model calculated with (a) WENO-5 and (b) 
WENO-5-AD. The temperature range is 5500 K (blue) to 11500 K (red). 
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Figure 18: A strong vortex tube reaching up to the high photosphere and, possibly there- 
fore, the chromosphere. The facing side is 1 Mm wide. Volume rendering: V(p — p). Slice: 
temperature. 
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